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Abstract 

We propose an extension of the nonequilibrium invaded cluster (IC) algorithm, which reestab- 
lishes a correct scaling of fluctuations at criticality and also self-adjusts to the critical temperature. 
We show that by introducing a single constraint to the intrinsic quantity of the IC algorithm the 
temperature becomes well defined and the sampling of the equilibrium ensemble is regained. The 
procedure is applied to the Potts model in two and three dimensions. 
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I. INTRODUCTION 



Monte Carlo (MC) studies of phase transitions have given rise to several cluster algorithms 
such as the Swendsen-Wang (SW) [l| or Wolff algorithm [2I, which brought an important 
improvement to simulations by reducing the critical slowing down. More recently another 
cluster algorithm was proposed by Machta et al. j^, inspired by the invasion percolation 
with an additional advantage that it self-regulates to the critical point and does not require 
prior knowledge of the critical temperature. It was applied to various classical models 
of phase transitions , and extended to more complex cases, from frustration [5| to 
tricritical points j^. The price to be paid is that the configurations generated by this self- 
adjusting nonequilibrium process are not distributed according to the canonical ensemble 







Thus, although the IC algorithm can be used to analyze a certain number of properties 
at criticality, those which are describing fluctuations remain out of its reach. This question 
was analyzed from different aspects {3, S] and an alternative self-adiusting algorithm based 
on a modification of the Swendsen-Wang algorithm was proposed 9|]. 

By following a different approach, we propose an equilibriumlike invaded cluster (EIC) 
algorithm obtained by constraining temperature uncertainty characteristic to the invaded 
cluster (IC) algorithm into limits compatible with the equilibrium distribution. We show 
that by applying this single constraint to the IC algorithm, correct scaling properties of 
thermodynamic observables are reproduced, while the algorithm still gives the critical tem- 
perature as an output. The method is demonstrated on the example of the Potts model in 
two and three dimensions. 

We consider the Potts model defined by the Hamiltonian 

H = -JJ2 i^s.,s,-l), (1) 

<i,j> 

where Sj denotes the g-state Potts variable at the lattice site i and the summation runs over 
the nearest neighbors only. The cluster algorithms including the IC algorithm are based 



on the Fortuin and Kasteleyn (FK) Uj expansion, which shows that the partition function 
of the model is equivalent to the one of the random-cluster (RC) model, which can be 
written as 

Z = ^/(^^l - pf~''^^k''^^\ (2) 
7er 
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and understood as a generalized bond percolation, where 

p=l- e-^' 



(3) 



is the bond probability {(3 is the Boltzmann factor). The summation in Eq. ([2]) runs over 
the set of all the graphs on the lattice F, while each graph 7 represents one possible bond 
configuration. B is the number of the lattice edges, 6(7) denotes the number of bonds, and 
0(7) is the number of connected components (FK clusters) in the graph 7. 



II. ALGORITHM 



In standard cluster algorithms such as the SW or Wolff algorithm, FK clusters are formed 
on a lattice by adding bonds between neighbors in the same state with a temperature- 
dependent probability p ([3]) . The configuration is then randomized by flipping all the clusters 
to different states and the procedure is repeated from the beginning. In the IC algorithm by 
Machta et al. the cluster formation phase is different since the probability p is not given in 
advance. Bonds are placed randomly between neighbors in the same state, until one of the 
clusters percolates. This is considered as a signature of the critical point in the system of 
finite size. The adding of bonds then stops, clusters are randomized and ready for the next 
iteration. The output that is recorded is the ratio 

where 6(7) is the number of added bonds and 6(7) the number of neighboring pairs in the 
same state for a given graph 7 on a lattice. The average p^ is identified with the quantity 
p, yielding a critical temperature estimate. The procedure is self-adjusting to the critical 
temperature imposed by the stopping condition. Namely, if during one step the excessive 
number of bonds is added, which would statistically contribute to a much lower temperature, 
in the next step the clusters will percolate easier due to large number of satisfied neighbors 
created, which will statistically contribute to temperatures much higher than critical. The 
problem is that the resulting oscillations in p^ remain too large during the entire process. 
As a result it does not sample the canonical ensemble, but its proper IC ensemble, giving 
much broader distribution. Since in the thermodynamic limit \Jvar{E) / E still holds, 
the IC algorithm is assumed to become equivalent in that limit to the canonical ensemble 
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and to produce the same average values, including the correct Tc- However, the convergence 
with size is more difficult to control. Also, the variances in energy or magnetization do not 
describe the equilibrium fluctuations, but rather the specific IC dynamics. 

It is important here to notice that the variable p of Eq. ([2]) is an a priori bond probability 
and consequently p^ should be distributed according to Gaussian distribution with p^ rc ~ P 



and the width \/var{p^)\Fic oc L~'^l'^ . In the equilibrium cluster algorithms p-y appears as a 
random input variable which obviously is distributed according to normal distribution. As 
already mentioned, the distribution of p^ generated as an output of the IC algorithm is far 
broader than L"^!'^ [J, I?]. 

We propose a simple restriction to the IC algorithm, which sets the allowed range of p^ 
to be proportional to L~'^l'^ and strictly less than \/ var{j)^)\Rc- The procedure is straight- 
forward and works as follows. The iterations are grouped in intervals of Na MC steps. 

(a) During the first interval of Na MC steps we follow the original IC rule. The average 
value p^ Q over the first Na iterations is found, to be used as an approximation for the next 
Na steps. 

(b) In the next interval of Na MC steps the stopping rule is modified by requiring the 
system to percolate but with the constraint p^ Q — v < p^ < p^^^ + v, where v is the parameter 
set to be proportional to L""^/^. If the system percolates before the lower bound of p-y is 
reached, the bonds are still added until the lower bound is attained. If the upper bound 
of p-y is reached before percolation, the process is stopped without requiring the system to 
percolate. At the end of the interval, the new average over p^, p^^-^ is calculated. 

(c) In every consecutive interval i of A''^ steps the bounds of p^ are set by — v < 
Py < py j_i + V and the same stopping rule as in (b) applied. 

The first few intervals of Na steps are discarded and the configurations are recorded after 
the steady state has been reached. 

We remark that the width of the resulting distribution of p^ has two contributions: (a) 
from the fiuctuations of Pj within a group of Na MC steps proportional to oc 
(b) from the fiuctuations of the mean value p^ ^ that correspond to the actual temperature 
fiuctuations, producing the width of order L~^/'^ < if the heat capacity exponent 

a > 0, or of order L""^/^ if a < 0. Even when it decays as L""^/^, it is much smaller than 
the contribution from (a) because is a result from averaging over Na MC steps. The 
value of Na is required to be large enough to allow the uncertainty of values of ^ to be less 
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than the variation of j between the groups i. The ^yvar(p^) of EIC algorithm thus can 
be regulated by changing v and is proportional to L~^/^. 

Consequently, the width of the distribution of allowed during the entire process of 
iterations is within the same limits as in the conventional cluster algorithms, and we can 
consider that the temperature during the whole process is equally well defined. 

The question of the ensemble that EIC algorithm generates is less trivial and we do 
not attempt here any detailed study of the exact connection between the EIC and the 
canonical ensemble. We only point out that there is a nonzero probability of generating 
any spin configuration with any value of p^, so the ergodicity of the procedure is fulfilled. 
The numerical evidence will be given below that the EIC algorithm samples the canonical 
ensemble well. 

III. RESULTS 

The EIC algorithm was applied to the cases g = 2, 3 and g = 2 of the Potts model in two 
and three dimensions, respectively. In two dimensions (2iD) exact results are known both 
for critical temperature and for the exponents, while in 3D very accurate estimations are 
available [l^ . The simulations were performed on lattices with periodic boundary conditions, 
and moderate sizes up to L = 226 and L = 64 in 2D and 3D, respectively. The statistics 
varied from 10^ iterations for smaller to 2 ■ 10^ iterations for the largest lattice sizes. Free 
parameters of the algorithm were chosen to be w = L~^^'^/10 and interval Na = 100 Monte 
Carlo steps (MCS). Percolation was established by the topological rule, i.e. by the condition 
that the infinite cluster wraps around the lattice. 

We include for comparison the results obtained by the SW algorithm with the same 
lattice sizes and using the same statistics performed at two different temperatures: (a) at 
infinite system critical temperature; i.e., for p = Pc{L oo); (b) at L-dependent critical 
temperatures; i.e. for pc{L) approximated by previously calculated Pc{L)\eic- Such a choice 
is justified because the difference between Pc{L)\eic and the position of e.g. the susceptibility 
maximum obtained by the SW algorithm for the lattice size L is much smaller than its width. 

In Fig. [Hare compared the energy histograms, generated by EIC, SW and IC algorithms 
for the case q = 2, D = 3. All three sets were taken for sizes L = 16 and L = 36 and rescaled 
with the same exponent x, set to produce the collapsing fit of the SW data. The distribution 
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Figure 1: (Color online) Energy histograms 'w{e) of the EIC, SW and IC algorithm for the 3D 
Ising model and sizes L = 16 and L = 36 rescaled with the same exponent x « 1.34. e denotes the 
averaged energy of each individual curve. Histograms of SW algorithm were taken at Pc{L)\eic- 



generated by EIC algorithm is the narrowest one, due to the choice of the arbitrary constant 
factor in the parameter v. More important is that the width of the EIC histogram scales 
with the same exponent as the one of the canonical ensemble obtained by SW algorithm, so 
that the fluctuations in both ensembles follow the same law. Also, by tuning the constant 
factor in v, very good overlap can be established between the EIC and SW histograms of all 
sizes. On the other hand, the width of the IC ensemble scales differently, consistent with the 
fact that the corresponding fluctuations follow the specific IC dynamics, and do not describe 
the equilibrium energy fluctuations at criticality. 

Let us present now the obtained results for critical point and critical exponents summa- 
rized in Table I. 

First we examine the convergence of Pc{L) with size. In Fig. [2] are presented data for 
Pc{L) vs L"^/*^, where the exact (or best known) value for was used. The obtained linear 
fits are very clean in comparison to the extrapolations within the IC algorithm (cf., e.g.. Fig. 
1 of the first Ref. in [sl), where the use of medians was necessary to reduce the additional 
size effects. The extrapolations to the limit L ^ oo are compared in Table I to the known 



results for pc in 2D (calculated from the analytic expression Pc = 1 — (1 + ^yq) ^ (l2|) and 



in 3D [l^ and show the precision up to four or five significant digits already for modest 
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sizes used here. The estimations of the temperature critical exponent, \/v were obtained 




Figure 2: (Color online) Data oipc{L) for cases I? = 2, g = 2, 3 and D = 3, g = 2. Lines represent 
linear fits, using exact or best known results for ^. 

from the power law form Pc{L) — pc{oo) oc L~^^'^ by taking the exact (best known) value for 
Pc(oo). 

The magnetic critical exponent was calculated in three different ways. We have calculated 



the fractal dimension of the percolation cluster s. 



max I p=Pc 



OC L^'^ and the anomalous dimension 



of the magnetization mL=p^ oc L^^'^ defined through the most populated of the q Potts states 



a. I.e. 



m 



We also examined the fluctuations of magnetization (see Fig [3]) which in equilibrium are 
related to the susceptibility x\p=pc L'^^" by the relation 



X 



(6) 



which assumes the validity of the fluctuation-dissipation theorem, not fulfilled for the stan- 
dard IC algorithm. 

EIC results for all the quantities presented in Table I agree with the exact, or best 
approximate results within a few percent. Although we omit here the analysis of convergence, 
visible improvement of results was observed with increasing L, so that discrepancies could 
be attributed to relatively modest sizes used. This is also supported by the fact that the 



Table I: Results for critical exponents obtained by EIC (from L = 64 to 226 for 2D and from 
L = 16 to 64 for 3D) compared to those by SW algorithm obtained for the same lattice sizes with 
the same statistics, and to the exact, or best approximate results. 





Pc Vt 


Vh 


V 


1 

V 


2D, g = 2 


EIC 


0.58575(2) 0.98(1) 


1.868(6) 


0.132(5) 


1.775(7) 


SW 




1.875(4) 


0.125(3) 


1.765(5) 


SW'' 




1.876(5) 


0.127(4) 


1.750(5) 


exact 


0.585786 ... 1 


15 

8 


1 

8 


7 
4 


2B,q = 3 


EIC 


0.63397(2) 1.19(2) 


1.861(6) 


0.140(5) 


1.745(7) 


SW 




1.86(5) 


0.136(5) 


1.750(8) 


SW'' 




1.88(4) 


0.120(4) 


1.740(5) 


exact 


0.633974 ... 1 


28 
15 


2 

15 


26 
15 


3D, g = 2 


EIC 


0.35809(1) 1.590(3) 


2.481(3) 


0.520(2) 


1.987(4) 


SW^ 




2.495(5) 


0.506(5) 


1.992(5) 


SW'' 




2.490(5) 


0.510(4) 


1.989(5) 


[13, i^"^ 


0.358098(3) 1.587(1) 


2.482(1) 


0.5181(5) 


1.963(1) 


"Simulations at Pc{L — s- oo). 
''Simulations a^^c{L)\Eic- 
'^Exact values [12]. 

''Best known values for critical point 13 1 and exponents 


14|. 



simulations by SW algorithm with the same sizes produced very similar results. Compared 
to those of the standard IC algorithm [sl, the fits of quantities derived from the mean values 
at criticality, such as pc, 1/^5 /^/^ significantly improved. Moreover, the correct 
exponent is regained for the fluctuations of the magnetization, and the susceptibility shows 
the correct scaling within the same error margins as within a SW algorithm. 
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ln(L) 



Figure 3: (Color online) A In-ln plot of the fluctuations of magnetization depending on the system 
size L. 

IV. DISCUSSION AND CONCLUSION 



To summarize, we have proposed an extension of the IC algorithm which recovers the 
correct sampling of the equilibrium ensemble at criticality and still preserves self-regulation 
to the critical temperature. By imposing on the distribution of the variable the width 
proportional to L~^^'^ (compatible with the Gaussian distribution), we reduced the uncer- 
tainty of the temperature variable to the one in the equilibrium algorithms and obtained 
the correct scaling of the fluctuations at criticality. Furthermore, our intervention does not 
slow down the IC algorithm, and the running times per MCS for the ETC algorithm remain 
the same. For example, a run of 10^ MCS for 2D Ising model on L = 64 lattice requires 
approximately 500 s on the AMD Opteron 240 processor (1.4 GHz). The procedure is il- 
lustrated on several cases of second-order phase transitions in the Potts model in two and 
three dimensions belonging to different universality classes. 

In comparison with the self- regulating algorithm by Tomita and Okabe , this approach 
is conceptually different: their approach is an extension of the SW algorithm which allows 
variation in p, which is directly related to the fluctuations of temperature, and the corre- 
sponding width is oc L~^/^ . As far as the efficiency is compared, the authors of Ref. 9| 
argue their algorithm to be faster per individual iterations since it does not require multiple 
checking of percolation during a single iteration. On the other hand, in Ref. 91] not less than 
10000 MCS preparation steps were required before the iterations could be recorded, while 
in the EIC approach not more than 2000 MCS were sufficient. Thus, it would be interesting 
to compare the autocorrelation times for the two methods. 



In the future, more detailed study remains to be done of the ensemble, the leading 
convergence exponent, the autocorrelations, and the corresponding dynamic exponent. 

The improved convergence and possibility to calculate the correct scaling of fluctuations 
may be useful in a number of problems. It may have particular advantage in the cases with 
quenched disorder and lack of self-averaging, where calculations at the sample-dependent 
critical temperatures have to be performed, and where the standard IC algorithm was of 



limited efficiency 
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